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Abstract 

We consider, by means of the variational approximation (VA) and direct numer- 
ical simulations of the Gross-Pitaevskii (GP) equation, the dynamics of 2D and 3D 
condensates with a scattering length containing constant and harmonically varying 
parts, which can be achieved with an ac magnetic field tuned to the Feshbach res- 
onance. For a rapid time modulation, we develop an approach based on the direct 
averaging of the GP equation, without using the VA. In the 2D case, both VA and 
direct simulations, as well as the averaging method, reveal the existence of stable 
self-confined condensates without an external trap, in agreement with qualitatively 
similar results recently reported for spatial solitons in nonlinear optics. In the 3D 
case, the VA again predicts the existence of a stable self-confined condensate with- 
out a trap. In this case, direct simulations demonstrate that the stability is limited 
in time, eventually switching into collapse, even though the constant part of the 
scattering length is positive (but not too large). Thus a spatially uniform ac mag- 
netic field, resonantly tuned to control the scattering length, may play the role of 
an effective trap confining the condensate, and sometimes causing its collapse. 
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1 Introduction 



Collisions between atoms play a crucially important role in the dynamics of Bose- 
Einstein condensates (BECs). As is commonly known, the collisions are accounted 
for by the cubic term in the corresponding Gross-Pitaevskii (GP) equation that 
describes BEG dynamics in the mean-field approximation. The coefficient in front 
of the cubic term, proportional to the collision scattering length, may be both 
positive and negative, which corresponds, respectively, to repulsive and attractive 
interactions between the atoms [1] . In the case of an attractive interaction, a soliton 
may be formed in an effectively one-dimensional (ID) condensate [2]; however, in 
2D and 3D cases the attraction results in the collapse of the condensate {weak and 
strong collapse, respectively [3]) if the number of atoms exceeds a critical value [1]. 

Recently developed experimental techniques [4] make it possible to effectively 
control the sign of the scattering length using an external magnetic field because 
the interaction constant can be changed through the Feshbach resonance [5]. This 
technique makes it possible to quickly reverse (in time) the sign of the interaction 
from repulsion to attraction, which gives rise, via the onset of collapse, to an abrupt 
shrinkage of the condensate, followed by a burst of emitted atoms and the formation 
of a stable residual condensate [4]. 

A natural generalization of this approach for controlling the strength and sign 
of the interaction between atoms and, thus, the coefficient in front of the cubic term 
in the corresponding GP equation, is the application of a magnetic field resonantly 
coupled to the atoms and consisting, in the general case of dc and ac components. 
The dynamical behavior of 2D and 3D condensates in this case is then an issue of 
straightforward physical interest, as it may be readily implemented in experiments. 
This is the subject of the present work. 

It is noteworthy that, in the 2D case, this issue is similar to a problem which was 
recently considered in nonlinear optics for (2-|-l)D spatial solitons (i.e., self-confined 
cylindrical light beams) propagating across a nonlinear bulk medium with a layered 
structure, so that the size [6] and, possibly, the sign [7] of the Kerr (nonlinear) 
coefficient are subject to a periodic variation along the propagation distance (it 
plays the role of the evolutional variable, instead of time, in the description of 
optical spatial solitons). The same optical model makes also sense in the (3+l)D 
case, because it applies to the propagation of "light bullets" (3D spatiotemporal 
solitons [8]) through the layered medium [7]. We will demonstrate below that the 
results obtained for the BEG dynamics in the GP equation involving both a dc and 
ac nonlinearity are indeed similar to findings reported in the framework of the above- 
mentioned optical model. To the best of our knowledge, a GP equation with a rapid 
time-periodic modulation of this type is proposed in this work for the first time. 
Previously, a quasi-lD model was considered in which the BEG stability was affected 
by a rapid temporal modulation applied to the trapping potential (rather than to 
the spatially uniform nonlinearity coefficient) [9] and the macroscopic quantum 
interference and resonances have been studied in [10]. Resonances in 2D and 3D 
BEG with periodically varying atomic scattering length has been considered in 
[11, 12]. 

The paper is organized as follows. In section 2, we formulate the model to be 
considered in this work and VA which will be employed to analyze the model. In 
section 3, variational and numerical results are presented for the 2D case (the anal- 
ysis based on VA also employs the Kapitsa averaging procedure). Both approaches 
demonstrate the existence of a stable self-sustained condensate, in a certain region 
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of parameter space, so that the condensate can be effectively confined and main- 
tained by means of a spatially uniform resonant ac magnetic field, without any 
trapping potential. In section 3, we also develop an alternative analytical approach, 
based on the application of the averaging procedure directly to the GP equation, 
without using the VA. Results produced by this approach confirm those obtained by 
means of VA. In section 4, we show that the results for the 3D case are essentially 
different from the ones in the 2D case. Here VA also predicts the possibility of a 
stable condensate, while direct simulations demonstrate that the stability is limited 
in time, finally giving way to collapse; a noteworthy fact is that, while VA per se 
still provides reasonable results in the 3D case, the averaging proc;edurc, if com- 
bined with VA, may yield completely wrong predictions in this case. A nontrivial 
feature demonstrated by direct simulations in the 3D case is that the ac component 
of the nonlincarity may give rise to collapse oven in the case when the dc (constant) 
component corresponds to repulsion. The paper is concluded by section 5. 



2 The model and variational approximation 

We take the mean-field GP equation for the single-particle wave function in its usual 
form, 

z^^ = -^AV. + ,H^^, (1) 

where A is the 2D or 3D Laplacian, r is the corresponding radial variable and 
g = Airh'^as/m, where as, m are respectively the atomic scattering length and mass. 
As indicated above we will assume the scattering length to be time-modulated so 
that the nonlinear coefficient in Eq. (1) takes the form g = go + gi sin(xi), where ao 
and ai are the amplitudes of dc and ac parts, and x is the ac-modulation frequency. 

Usually an external trapping potential is included to stabilize the condensate. 

We have omitted it because it does not play an essential role. This is also the case in 
some other situations, e.g., the formation of a stable Skyrmion in a two-component 
condensate [13]. In fact, we will demonstrate that the temporal modulation of the 
nonlinear coefficient, combining the dc and ac parts as in Eq. (2) may, in a certain 
sense, replace the trapping potential. Another caveat concerning the present model 
is that the frequency of the ac drive must be chosen far enough from resonance 
with any transition between the ground state of the condensate and excited states, 
otherwise the mean-field description based on the GP equation will not be adequate. 

We now cast Eq. (1) in a normalized form by introducing a typical frequency 
SI ~ 2gno/fi, where no is the condensate density and rescale the time and space 
variables as t' = fit r' = r-^/2mfl/h. This leads to the following equation where the 
' have been omitted 

in which it is implied that depends only on t and r, £> = 2 or 3 is the spatial 
dimension, Ao,i = —go,i/ {S^K), a; = x/fi. 

Note that Ao > and Ao < in Eq. (2) correspond, respectively, to the self- 
focusing and self-defocusing nonlinearity. Rescaling the field tp, we will set |Ao| = 1, 
so that Ao remains a sign-defining parameter. 
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The next step is to apply the VA to Eq. (2). This approximation was originally 
proposed [14] and developed in nonlinear optics, first for ID problems and, later 
for multi-dimensional models (see a recent review [15]). A similar technique was 
elaborated for the description of the multidimensional BEC dynamics based on the 
GP equation. [16]. 

To apply VA in the present case, we notice that the Lagrangian density gener- 
ating Eq. (2) is 



)- 






dr 



+ ^A(^)|V|^ (3) 



where X{t) = Ao + Ai sm{u)t), and the asterisk stands for the complex conjugation. 
The variational ansatz for the wave function of the condensate is chosen as the 
Gaussian [14], 

Mr,t)=Ait)exp(^-^^ + ^tb{t)r^+iS{t)^ , (4) 

where A, a, b and S are, respectively, the amplitude, width, chirp and overall phase, 
which are assumed to be real functions of time. We did not include the degree 
of freedom related to the coordinate of the condensate's center, as the trapping 
potential, although not explicitly included into the model, is assumed to prevent 
the motion of the condensate as a whole. 

Following the standard procedure [15], we insert the ansatz into the density (3) 
and calculate the effective Lagrangian, 



/■OO 

= Cd / C{i;g)r''-'dr, 
Jo 



(5) 



where Cd = 27r or An in the 2D or 3D cases, respectively. Finally the evolution 
equations for the time-dependent parameters of the ansatz (4) are derived from L^s 
using the corresponding Eulcr - Lagrange equations. Subsequent analysis, as well 
as the results of direct numerical simulations, are presented separately for the 2D 
and 3D cases in the two following sections. 



3 The two-dimensional case 
3.1 Variational approximation 

In the 2D case, the calculation of the effective Lagrangian (5 ) yields 

Lf/) = n{~a^A% - a^AH - A^ - a^AH^ + ^X{t) a^A^), (6) 

where the overdot stands for the time derivative. The Euler-Lagrange equations 
following from this Lagrangian yield, the conservation of the number of atoms A'' in 
the condensate, 

ttA^o^ = N = const, (7) 
an expressions for the chirp and the width. 




X{t)N 
27ro4 ' 
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and a closed-form evolution equation for the width: 

rf^a 2(2 - X{t)N/2n) 



which we rewrite as 



(fa — A + esin(a;t) 

dt2" = ^3 ' 



(8) 



(9) 



A = 2 (AoA^/(27r) - 2) , e = -XiN/w. (10) 

In the absence of an ac component, e = 0, Eq. (9) conserves the energy i?2D = 
(d^ — Aa~^) /2. Obviously, E2U —00 as a — » 0, if A > 0, and i?2D +00 
as a — > 0, if A < 0. This means that, in the absence of the ac component, the 
2D pulse is expected to collapse if A > 0, and to spread out if A < 0. The case 
A = corresponds to the critical number of particles in the condensate (the so-called 
"Townes soliton"). Note that a numerically exact value of the critical number is 
(in the present notation) N = 1.862 [3], while the variational equation (10) yields 
= 2 (if Ao = +1). 

It is natural to specially consider the case when the ac component of the non- 
linear coefficient oscillates at a high frequency. In this case, Eq. (9) can be 
treated analytically by means of the Kapitsa averaging method. To this end, we set 
a{t) = a+da, with \da\ « |a|, where a varies on a slow time scale and Sa is a rapidly 
varying function with a zero mean value. After straightforward manipulations, we 
derive the following equations for the slow and rapid variables, 

-^a = -A(a-3 + 6a-5 ((Ja^)) - 3e {6a sin(w<)) a'^, (11) 

-^5a = 3 5aAd"^-Fesin(a;f)d"^. (12) 

where (...) stands for averaging over the period 277/0;. A solution to Eq. (12) is 

r / X esin(wi) 

the substitution of which into Eq. (11) yields the final evolution equation for the 
slow variable, 



d2 3 



-A 



3Ae2 



(a;2a4 + 3A)2 2L0^d^ + 3A 



(14) 



To examine whether collapse is enforced or inhibited by the ac component of 
the nonlinearity, one may consider Eq. (14) in the limit d — > 0. In this limit, the 
equation reduces to 

^aM-A + ii)S- (15) 

It immediately follows from Eq. (15) that, if the amplitude of the high-frequency ac 
component is large enough, > 6A^ , the behavior of the condensate (in the limit of 
small d) is exactly opposite to that which would be expected in the presence of the 
dc component only: in the case A > 0, bounce should occur rather than collapse, 
and vice versa in the case A < 0. 

On the other hand, in the limit of large d, Eq. (14) takes the asymptotic 
form d^d/dt^ = — Ad~^, which shows that the condensate remains self-confined 
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in the case A > i.e., if tlic number of atoms exceeds the critical value. This 
consideration is relevant if a though being large remains smaller than the limit 
imposed by an external trapping potential, should it be added to the model. Thus, 
these asymptotic results guarantee that Eq. (14) gives rise to a stable behavior of 
the condensate, both the collapse and decay (spreading out) being ruled out if 

e>v^A>0. (16) 

In the experiments with for example with ^Li with the critical number ~ 1500 
atoms if we have initially 1800 atoms (i.e. N/2'k = 2.2) to stabilize the condensate 
this means that the atomic scattering length for Aq = 1 should be harmonically 
modulated with the amplitude e = 0.98. In fact, the conditions (16) ensure that 
the right-hand side of Eq. (14) is positive for small a and negative for large a. This 
implies that Eq. (14) must give rise to a stable fixed point (FP). Indeed, when the 
conditions (16) hold, the right-hand side of Eq. (14) vanishes at exactly one FP, 

which can be easily checked to be stable through the calculation of an eigenfrequency 
of small oscillations around it. 

Direct numerical simulations of Eq. (9) produce results (not shown here) which 
are in exact correspondence with those provided by the averaging method, i.e., a 
stable state with a{t) performing small oscillations around the point (17). The 3D 
situation shows a drastic difference because this correspondence breaks down, as 
shown in the next section. 

For the sake of comparison with the results obtained by means of an alternative 

approach in the next subsection, wc also need an approximate form of Eq. (14) 
valid in the limit of small A (i.e., when the number of atoms in the condensate is 
close to the critical value) and very large ui: 

To estimate the value of the amplitude of the high-frequency ac component 

necessary to stop the collapse, we note that a characteristic trap frequency is ~ 
100 Hz. So, for a modulation frequency ~ 3 kHz, which may be regarded as a 
typical "high modulation frequency", the dimensionless w is ~ 30. If the initial 
dimensionlcss number of atoms is, for example, N/2tt = 2.2 so that according to 
Eq. (10), A = 0.4 (this corresponds to the ^Li condensate with « 1800 atoms, the 
critical number being « 1500), and the parameters of modulation arc Aq = 1, Ai = 
2.3, e = 10, then the stationary value of the condensate width found from (17) is 
ttst = 0.8Z, where I = y^mil/h is the healing length. 

Thus our analytical approach, based on the VA and the subsequent use of the 
assumption that the number of atoms slightly exceeds the critical value, leads to an 
important prediction: in the 2D case, the ac component of the nonlincarity, acting 
jointly with the dc one corresponding to attraction, may give rise not to collapse, 
but rather to a stable soliton-like oscillatory condensate state which confines itself 
without the trapping potential. It is relevant to mention that a qualitatively similar 
result, viz., the existence of stable periodically oscillating spatial cylindrical solitons 
in a bulk nonlinear-optical medium consisting of alternating layers with opposite 
signs of the Kerr coefficient, was reported in Rcf. [7], where this result was obtained 
in a completely analytical form on the basis of the VA, and was confirmed by direct 
numerical simulations. 
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3.2 Averaging of the Gross-Pitaevskii equation and Hamil- 
tonian 



In the case of a high-frequency modulation, there is a possibihty to apply the av- 
eraging method directly to the 2D equation (2), without using the VA. Note that 
direct averaging was applied to the 2D nonlinear Schrodinger equation (NLS) with 
a potential rapidly varying in space, rather than in time, in Ref. [17], where the 
main results were a renormalization of the parameters of the 2D NLS equation and 
a shift of the collapse threshold. As we will see below, a rapid temporal modulation 
of the nonlinear term in the GP equation leads to new effects, which do not reduce 
to a renormalization. Namely, new nonlinear-dispersive and higher-order nonlinear 
terms will appear in the corresponding effective NLS equation [see Eq. (22) below] . 
These terms essentially affect the dynamics of the collapsing condensate. 

Assuming that the ac frequency to is large, we rewrite the 2D equation (2) in a 
more general form, 

idxp/dt + A'ilj + X{ujt)\ip\'^ilj = 0, (19) 

where A is the 2D Laplacian. To derive an equation governing the slow variations 
of the field, we use the multiscale approach, writing the solution as an expansion in 
powers of 1 /uj and introducing the slow temporal variables, Tk = w~''t, fc = 0, 1, 2... , 
while the fast time is C = wt. Thus, the solution is sought for as 

V(r,i) = A{r,Tk)+oj-\i{C,A)+oj-^U2{C,A) + (20) 

with (uk) — 0, where (...) stands for the average over the period of the rapid mod- 
ulation, and we assume that Ao = -|-1 (i.e., the dc part of the nonlinear coefficient 
corresponds to attraction between the atoms). 

Following a procedure developed, for a similar problem, in Ref. [18], we first 
find the first and second corrections, 

ui = - (mi)] 1^1' A Ml = t [A(t) - (Ai)] dr, (21) 

Jo 



U2 = [M2 



{^i2)mA\^At + zA^A; 



A{\A\^A)] 



\A\'A[-[{f,, - (Ml))' - 2M] + (A) (m2 - (m2))]. 

HereM2 = /o^(mi- < Mi >)ds,M = (l/2)(< m? > - < Mi >') = (l/2)(< A^ > -1) 
(recall we have set |Ao| = 1). Using these results, we obtain the following evolution 
equation for the slowly varying field A{x,To), derived at the order oj~^: 

dA 



i^ + AA + \AfA + 2M(-) [\AfA- 



3\A\'^AA + 2\A\'^Ai\A\'^A) + A^Al"^ A*] = 0, (22) 

where e is the same amplitude of the ac component as in Eq. (10). We stress that 
Eq. (22) is valid in both 2D and 3D cases. In either case, it can be represented in 
the quasi-Hamiltonian form. 



i+6M(i)^Ar 



dA 
'dt 



.5Hq 

''SA*' 



H, = jdV [|VA|' - 2M (1)' \A\' - \ \At + 4M (1)' |v (|A|' A) 



(23) 



(24) 
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where dV is the infinitesimal volume in the 2D or 3D space. To cast this result 
in a canonical Hamiltonian representation, one needs to properly define the cor- 
responding symplectic structure (Poisson's brackets), which is not our aim here. 
However, we notice that, as it immediately follows from Eq. (23) and the reality of 
the (quasi-)Hamiltonian (24), Hq is an integral of motion, i.e., dHq/dt = 0. 

For a further analysis of the 2D case, we apply a modulation theory developed 
in Ref. [19]. According to this theory, the solution is searched for in the form of 
a modulated Townes soliton. The (above-mentioned) Townes soliton is a solution 
to the 2D NLS equation in the form ip{r,t) = e^^Rxir), where the function Rt{t) 
satisfies the boundary value problem 



Rj- 



^ R'j- — Rt + Rj- 



0, RIt{0) = 0, Rt{oo) = 0. 



(25) 



For this solution, the norm N and the Hamiltonian H take the well-known values, 

^2 1 



pOO pC 

Nt= R!^{r)rdr = AT^ = 1.862, Ht = 
Jo Jo 



{R'tY - ^Rrir) 



rdr = 0. (26) 



The averaged variational equation (22) indicates an increase of the critical num- 
ber of atoms for the collapse, as opposed to the classical value (26). Using the 
relation (20), we find 



/■OO 

Ncrit= / IVl'rrfr = TVt + 2M(-)%, 
Jo ^ 



where Iq = 11. 178 (see Appendix 1). This increase in the critical number of atoms 
is similar to the well-known energy enhancement of dispersion-managed solitons in 
optical fibers with periodically modulated dispersion [20, 21]. 

Another nontrivial perturbative effect is the appearance of a nonzero value of 
the phase chirp inside the stationary soliton. We define the mean value of the chirp 
as 

J,-'lm(^r)rdr 
b = . 

Making use of the expression (21) for the first correction, we find 



X°^rdrRHRr-m5)J^'-'^drR'^ 
w ^' ^ J+°°r^drR^ 



b = —BM (/.I - (Ml)) , B = 3^ ' ' ^\ = 0.596. 



To develop a general analysis, we assume that the solution with the number 
of atoms close to the critical value may be approximated as a modulated Townes 
soliton, i.e. 

. 2 

A{r, t) « [a{t)]-^ RT{r/ait))e'^, S = a{t) + ^,a = a-^ (27) 

with some function a{t) (where the overdot stands for d/dt). If the initial power 
is close to the critical value, i.e., when jA' — A^d << Nc and the perturbation is 
conservative,!. e. 

Im j dV[A*F{A)] = 

as in our case, a method worked out in Ref. [19] makes it possible to derive an 
evolution equation for the function a{t), starting from the approximation (27). The 
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equation of modulation theory for width is 



"'"" = -^° + 4A%;;^^^('^' ^^^^ 

where 

^° - " " Mo ' 

and Mo = (1/4) r^drR\. w 0.55. The auxiliary function is given by 

h{t) = 2a(t)Re[^ ^ dxdyF(AT)e-^^(ii!T + pVRt{p))]. (29) 

In the lowest-order approximation, the equation takes the form (for the harmonic 
modulation) 

$ = -^ + ^- (30) 
where Ai = (TV - 7Vc)/Mo - Ce"^ /{uj'^at) and C is 

Values of integrals are given in the Appendix 1. Thus the averaged equation predicts 
the arrest of collapse by the rapid modulations of the nonlinear term in the 2D GP 
equation. The comparison of Eq. (30) with its counterpart (18), which was derived 
by means of averaging the VA-generated equation (9), shows that both approaches 
lead to the same behavior near the collapse threshold. The numerical coefficients 
in the second terms are different due to the different profiles of the Gaussian and 
Townes soliton. In this connection, it is relevant to mention a recent work [22], 
which has demonstrated that, generally, one may indeed expect good agreement 
between results for 2D solitons produced by VA and by the method based on the 
modulated Townes soliton. 

Let us estimate the value of the fixed point for the numerical simulations per- 
formed in Ref. [6] . In this work the stable propagation of soliton has been observed 
for two step modulation of the nonlinear coefficient in 2D NLSE. The modula- 
tion of the nonlinear coefficient was A = l-|-eifT>i>0, and A = 1 — e for 
2T > t > T. The parameters in the numerical simulations has been taken as 
T = e = 0.1, A^/(27r) = 11.726/(27r), with the critical number as = 11.68/(27r). 
The map strength M is M = e^T^/24. For this values we have Cc = 0.49, that 
agreed with the value ac ~ 0.56 following from the numerical experiment. 

Instead of averaging Eq. (2), one can apply the averaging procedure, also based 
on the representation (20) for the wave function, directly to the Hamiltonian of Eq. 
(2). As a result, the averaged Hamiltonian is found in the form 

H= f dxdyllVAl"" + 2M(-)2|V(|A|2A)|2 - hA\^ - 6M(-)2|A|8]. (32) 

A possibility to stop the collapse, in the presence of a rapid periodic modulation 
of the atomic scattering length, can be explained on the basis of this Hamiltonian. 
To this end, following the pattern of the usual virial estimates [3], we note that, 
if a given field configuration has compressed itself to a spot with a size p, where 
the amplitude of the A-field is ~ K, the conservation of the number of particles N 
[which may be applied to the A-field through the relation (20)] yields the relation 

~ N (33) 
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(recall D is the space dimension). On the other hand, the same estimate for the 
strongest collapse-driving and collapse-stopping terms [the fourth and second terms, 
respectively, in expression (32)] i?_ and H+ in the Hamiltonian yields 

(34) 

Eliminating the amplitude from Eqs. (34) by means of the relation (33), we conclude 
that, in the case of the catastrophic self-compression of the field in the 2D space, 
p — > 0, both terms take the same asymptotic form, , hence the collapse 
may be stopped, depending on details of the configuration. However, in the 3D 
case the collapse-driving term diverges as while the collapse-stopping term has 
the asymptotic form p~^, for p — > 0, hence in this case the collapse, generally 
speaking, cannot be prevented. 

Lastly, it is relevant to mention that, although the quasi-Hamiltonian (24) is not 
identical to the averaged Hamiltonian (32), the virial estimate applied to Hq yields 
exactly the same result: the collapse can be stopped in the 2D but not in the 3D 
situation. 



3.3 Direct numerical results 

The existence of stable self-confined soliton-like oscillating condensate states, pre- 
dicted above by means of analytical approximations for the case (16), when the 
dc part of the nonlinearity corresponds to attraction between the atoms, and the 
amplitude of the ac component is not too small, must be checked against direct sim- 
ulations of the 2D equation (2). In fact, it was quite easy to confirm this prediction 
[in the case Aq = — 1, i.e., when the dc component of the nonlinearity corresponds 
to repulsion, the direct simulations always show a decay (spreading out) of the 
condensate, which also agrees with the above predictions]. 

A typical example of the formation of a self-confined condensate, supported by 
the combination of the self-focusing dc and sufficiently strong ac components of the 
nonlinearity in the absence of an external trap, is displayed in Fig. 1. On the left 
panel we show the pulse collapse at t « 0.3 in the absence of modulation. In the 
presence of modulation the pulse is stabilized for about 40 periods after which it 
decays. Note the presence of radiation as the pulse adjusts to the modulation. 



4 The three-dimensional case 

4.1 The variational approximation and averaging 

The calculation of the effective Lagrangian (5) in the 3D case yields 

(35) 



-^^a^-2^+^A(.)^^-J-3.V 



cf. Eq. (6). The Euler-Lagrangc equations applied to this Lagrangian yield the 
mass conservation, 

^3/2^2^3 = ^ const, 

an expressions for the chirp, 

d = 2a6, 6=4-26^- '^^^^ 



2x/27r3/2a5 ' 
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and the evolution equation for the width of the condensate, 

(fa _ 4 X{t) N 

Note the difference of Eq. (36) from its 2D counterpart (8). 



(36) 



As in the 2D case, we renormalize the amphtiiclcs of the do and ac components 
of the nonlinearity, A = 2-^/'^Tr-^/^XoN and e = -2-^/'^Tr-^/'^XiN, and cast Eq. 
(36) in the normaUzed form, 

(Pa _ 4 -A + esin(a;t) 

In the absence of the ac term, e = 0, Eq. (37) conserves the energy 

E3D = + 2o-2 - ^Aa-3. 

Obviously, £^3d — oo as a ^ 0, if A > 0, and E-^q +oo if A < 0, hence one will 
have collapse or decay (spreading out) of the pulse, respectively, in these two cases. 

Prior to applying the averaging procedure (as it was done above in the 2D 
case), we solved Eq. (37) numerically, without averaging, to show that (within the 
framework of VA) there is a region in parameter space where the condensate, that 
would decay under the action of the repulsive dc nonlinearity (A < 0), may be 
stabilized by the ac component of the nonlinearity, provided that its amplitude is 
sufficiently large. To this end, we employed the variable-step ordinary differential 
equation (ODE) solver D0PRI5 [23], which is a combination of the Runge-Kutta 
algorithm of the fourth and fifth orders, so that the instantaneous truncation error 
can be controlled. 

In Fig. 2 we show the dynamical behavior of solutions to Eq. (37), in terms of 
the Poincare section in the plane (a, d), obtained for A — —1, e = 100, uj ~ lO"' • tt, 
and initial conditions a(t = 0) = 0.3, 0.2, or 0.13 and a(t = 0) = 0. As it is obvious 
from Fig. 2, in all these cases the solution remains bounded and the condensate 
does not collapse or decay, its width performing quasi-periodic oscillations. 

In fact, the corresponding stability region in the parameter plane (w/tt, e) is 
small, see Fig. 3. It is also seen from Fig. 3 that the frequency and amplitude of 
the ac component need to be large to yield this stability. Notice that, for frequencies 
larger than 10^ • tt, the width of the condensate a(t) assumes very small values in the 
course of the evolution (as predicted by VA) so that collapse may occur in practice 
for the solution of the full equation (2). 

The stability is predicted by VA only for A < 0, i.e., for a repulsive dc component 
of the nonlinearity. In the opposite case, the VA predict solely collapse. 

As Lo is large enough in the stability region shown in Fig. 3, it seems natural to 
apply the Kapitsa's averaging method to this case too. Doing it the same way as 
was described in detail in the previous section for the 2D case, we find the rapidly 
oscillating correction (5a(t), cf. Eq. (13), 

esin(a;t)d 

w2a5-12a + 4A ^ ' 

and then arrive at the evolution equation for the slow variable d(t) [cf. Eq. (14)]: 



da 4 



4a - A -h ^TTii — . -h e 



2e^ 2 6d-5A 



w^oS - 12d -I- 4A (a;2a5 _ + 4A)2 



(39) 
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In the limit a — > 0, Eq. (39) takes the form 



fa 




(40) 



cf. Eq. (15). Equation (40) predicts one property of the 3D model correctly, 
viz., in the case A < and with a sufficiently large amplitude of the ac component 
[e > (4/-\/3) |A|, as it follows from Eq. (40], collapse takes place instead of spreading 
out. However, other results following from the averaged equation (39) are wrong, 
as compared to those following from the direct simulations of the full variational 
equation (37), which were displayed above, see Figs. 2 and 3. In particular, a 
detailed analysis of the right-hand side of Eq. (39) shows that it does not predict 
a stable FP for A < 0, and docs predict it for A > 0, exactly opposite to what was 
revealed by the direct simulations. This failure of the averaging approach (in stark 
contrast with the 2D case) may be explained by the existence of singular points 
in Eqs. (38) and (39) (for both A > and A < 0), at which the denominator 
w^a^ — 12a + 4A vanishes. Note that, in the 2D case with A > 0, for which the 
stable state was found in the previous section [see Eq. (16)], the corresponding 
equation (14) did not have singularities. 

4.2 Direct simulations of the Gross-Pitaevskii equation in 
the three-dimensional case 

Verification of the above results given by VA against direct simulations of the 3D 
version of the radial equation (2) is necessary. The partial differential equation 
simulations were carried out by means of the method of lines implemented with the 
D0PRI5 ODE solver and space discretization involving high order finite differences, 
see the details in Appendix 2. The relative error in the conservation of the number 
of atoms was limited by 10~^. In the absence of the ac modulation, the energy was 
conserved with a relative error limited by 10~^. 

Quite naturally, in the case e = (no ac component) and A < 0, the simulations 

show straightforward decay of the condensate (not displayed here, as the picture 
is rather trivial). If an ac component of sufficiently large amplitude is added, sta- 
bilization of the condensate takes place temporarily, roughly the same way as is 
predicted by the solution of the variational equation (37). However, the stabiliza- 
tion is not permanent: the condensate begins to develop small-amplitude short-scale 
modulations around its center, and after about 50 periods of the ac modulation, it 
collapses. 

An example of this behavior is displayed in Fig. 4, for which iV = 1, A = — 1 
and io = 10** • tt. Figure 4 shows radial profiles of the density |M(r)p at different 
instants of time. 

Results presented in Fig. 4 turn out to be quite typical for the 3D case with 
A < 0. The eventual collapse which takes place in this case is a nontrivial feature, 
as it occurs despite the fact that the dc part of the nonlinearity drives the con- 
densate towards spreading out. Therefore, a basic characteristic of the system is 
a dependence of the minimum ac amplitude e, which gives rise to the collapse at 
fixed A = — 1, versus the ac frequency w. Several points marked by stars show this 
dependence in Fig. 3. It is quite natural that the minimum value of e necessary for 
the collapse grows with w. On the other hand, for oj not too large, the minimum 
ac amplitude necessary for the onset of collapse becomes small, as even a small e is 
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sufficient to push the condensate into collapse during the relatively long half-period 
when the sign of the net nonlinearity coefficient X{t) is positive, see Eq. (19). 

In the case of A > we have never been able to prevent the collapse of the pulse. 
This is in agreement with the analysis developed in the previous section on the basis 
of the Hamiltonian of the averaged version of the GP equation, which showed that 
the collapse cannot be stopped in the 3D case, provided that the amplitude of the ac 
component is large enough. Besides that, this eventual result is also in accordance 
with findings of direct simulations of the propagation of localized 3D pulses in the 
above-mentioned model of the nonlinear optical medium consisting of alternating 
layers with opposite signs of the Kerr coefficient: on the contrary to the stable 2D 
spatial solitons [7], the 3D spatiotemporal "light bullets" can never be stable in this 
model [24]. 



5 Conclusion 

In this work, we have studied the dynamics of 2D and 3D Bose-Einstein conden- 
sates in the case when the scattering length in the Gross-Pitaevskii (GP) equation 
contains constant (dc) and time-variable (ac) parts. This may be achieved in the 
experiment by means of a resonantly tuned ac magnetic field. Using the varia- 
tional approximation (VA), simulating the GP equation directly, and applying the 
averaging procedure to the GP equation without the use of VA, we have demon- 
strated that, in the 2D case, the ac component of the nonlinearity makes it possible 
to maintain the condensate in a stable self-confined state without external traps, 
which qualitatively agrees with recent results reported for spatial solitons in non- 
linear optics. In the 3D case, VA also predicts a stable self-confined state of the 
condensate without a trap, provided that the constant part of the nonlinearity cor- 
responds to repulsion between atoms. Direct simulations reveal that, in this case, 
the stability of the self-confined condensate is limited in time. Eventually, collapse 
takes place, despite the fact that the dc component of the nonlinearity is repulsive. 
Thus, we conclude that the spatially uniform ac magnetic field, resonantly tuned 
to affect the scattering length, may readily play the role of an effective trap which 
confines the condensate, and sometimes enforces its collapse. These predictions can 
be verified in experiments. 
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6 Appendix 1: calculation of the moments 

For the modulation analysis of section 3.2, we introduce the following integrals 
involving the Townes soliton. 

The boundary value problem (25) has been solved by discretzing using finite 
differences and using the shooting method. The solutions give a residual smaller 
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than 10"''. 



The integrals have been calculated using the trapezoidal rule. As a test the 
following integrals has been calculated, the norm N{Rs) = 1.862... and hamiltonian 
{H{Rg) = 0). For the other integrals we obtain 



h = I r'^drR^ = 1-7 ; l2= rdrR^{R'j.f = 2.529 
Jo Jo 



oo 



J3 = / rdrR^iR'TY = 5.730; I4, = / r^drR^R'Ty = -3.109 ; 
Jo Jo 

pOO pOO POO 

h = i drR^ = 11.472 ; 4 = / rdrR^ = 11.312 ; /g = / rdrR^ = 39.963, 
Jo Jo Jo 

POO POO 

Ig= rdrR^{R;rf = -4.872 ; ho = rdrR^{R'rpf = 3.669 ; 
Jo Jo 

POO 

hi = / rdrR^{Rirf = -2.314. 
Jo 

7 Appendix 2: numerical procedure for solving 
the partial differential equation 

Following [19], we have solved the cylindrical NLS equation (2) using the method of 

lines where the solution is advanced in time using an ordinary differential equation 
(ODE) solver and the spatial part is discretized using finite differences. Because of 
its implicit character, this method allows for great stability and accuracy as well 
as giving the possibility of implementing directly the cylindrical Laplacian and its 

associated boundary conditions. 

Specifically we use as ODE solver the variable step Runge-Kutta of order 4-5 
D0PRI5 [23]. which enables to control the error made at each step and bound it 
by a given tolerance. For all the runs presented the relative error at each step is 
below 10~^. The cylindrical Laplacian + {D — l)dr/r is approximated at each 
node n of the grid using the following formulas 

where h is the mesh size. We have therefore a method to solve (2) that is 0{dt'^, K^) 

The first node corresponds to r = and to its left we introduce two fictitious 
points so that ^l^r-r = = 0. At the right hand side boundary chosen sufficiently 
far from the pulse, V was set to be 0, again in two points. 

The number of mesh points was 4000 and the tolerance of the integrator set to 

10~^. In all cases the norm was conserved up to 10~* in relative error as 
was the Hamiltonian in the absence of modulation. The latter quantity provided 
an accurate indicator of collapse. 



14 



References 



[1] F. Dalfovo, St. Giorgini, L.P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 
463 (1999). 

[2] L. Khaykovich, F. Schrek, F. Ferrari, T. Bourdel, J. CubizoUes, L.D. Carr, 
Y. Castin, and C. Salomon, Science 296, 1290 (2002); K.E. Strecker, G.B. 
Partrage, A.G. Truscott, and R.G. Hulet, Nature 417, 150 (2002). 

[3] L. Bcrgc, Phys. Rep. 303, 259 (1998). 

[4] S.L. Cornish, N.R. Claussen, J.L. Roberts, E.A. Cornell, and C.E. Wieman, 
Phys. Rev. Lett. 85, 1795 (2000). 

[5] Y. Kagan, E.L. Surkov, and G.V. Shlyapnikov, Phys. Rev. Lett. 79, 2604, 1997. 

[6] L. Berge, V.K. Mezentsev, J. Juul Rasmussen, P.L. Christiansen, and Yu.B. 
Gaididei, Opt. Lett. 25, 1037 (2000). 

[7] I. Towers and B.A. Malomed, J. Opt. Soc. Am. 19, 537 (2002). 

[8] Y. Silberberg, Opt. Lett. 15, 1281 (1990); B.A. Malomed, P. Drummond, H. 
He, A. Berntson, D. Anderson, and M. Lisak, Phys. Rev. E 56, 4725 (1997). 

[9] R. Dum, A. Sampcra, K. Suominen, M. Brewczyk, M. Kus, K. Rzazevski, and 
M. Lewenstein, Phys. Rev. Lett. 90, 3899 (1998). 

[10] F.Kh. AbduUaev and R.A. Kraenkel, Phys.Lett. A 272, 395 (2000). 

[11] F.Kh. AbduUaev, J. Bronski, and R. Galimzyanov, eprint:cond-mat/0205464 
(2002). 

[12] S. Adhikari, cprint: cond-mat/0207579 (2002). 

[13] R.A. Battye, N.R. Cooper, and P.M. Sutcliffe, Phys. Rev. Lett. 88, 080401 
(2002). 

[14] D. Anderson, Phys. Rev. A 27, 1393 (1983). 
[15] B.A. Malomed, Progress in Optics 43, 69 (2002). 

[16] V.M. Perez-Garcia, H. Michinel, J.I. Cirac, M. Lewenstein, and P. ZoUer, Phys. 
Rev. A 56, 1424 (1997). 

[17] Yu. Kivshar, S. Turitsyn, Phys. Rev. E 49 R2536 (1994). 

[18] T.S. Yang and W.L. Kath, Opt. Lett. 22, 985 (1997). 

[19] G. Fibich and G.C. Papanicolaou, Phys. Lett. A 239, 167 (1998); SIAM Jour- 
nal of Appl. Math. 60, 183 (1999). 

[20] J.H.B. Nijhof, N.J. Doran, W. Forysiak and F.M. Knox, Electron. Lett. 33, 

1726 (1997). 

[21] F. Kh. AbduUaev and J.G. Caputo, Phys. Rev. E, 58, 6637, (1998). 

[22] L. Berge and A. Couairon, Physica D 152-153, 752 (2001). 

[23] E. Haircr, S. P. Norsctt and G. Wanner. Solving ordinary differential equations 
/ (Springer- Vcrlag, 1987). 

[24] I. Towers and B.A. Malomed, unpublished. 



15 



Figure 1: A typical example of the formation of a self-confined condensate, revealed 
from direct simulations of Eq. (2) in the two-dimensional case. The left panel shows 
pulse collapse in the absence of modulation for t ~ 0.3. The right panel shows the 
modulated pulse with the same initial condition for t « 0.6. The parameters are 
Ao = 2.4, Ai = 0.85 w = IOOtt and N = 5. 



Figure 2: The Poincarc section in the plane (a, a) for A = — l,e = 100, w = IO^tt, 
generated by the numerical solution of the variational equation (37) with different 
initial conditions (see the text). 



Figure 3: The region in the (e, w/tt) parameter plane where the numerical solution of 
Eq. (37) with A = — 1 predicts stable quasipcriodic solutions in the 3D case. Crosses 
mark points where stable solutions were actually obtained. Stars correspond to the 
minimum values of the ac-component's amplitude e eventually leading to collapse 
of the solution of the full partial differential equation (2) with A = — 1, see below. 



Figure 4: Time evolution of the condensate's shape |M|^(r) in the presence of the 
strong and fast ac modulation (w = IO'^tt, e = 90). From left to right the profiles of 
u^(r) are shown at times t = 0.007, 0.01 and 0.015. 



16 



6 





0.1 



0.2 



0.3 



0.4 



1000 



100 - 



10 



1 




10 



1000 



100000 



1e+07 



Omega/pi 



